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We discuss the use of a Langevin equation with a colored (correlated) noise to perform constant- 
temperature molecular dynamics simulations. Since the equations of motion are linear in nature, it 
is easy to predict the response of a Hamiltonian system to such a thermostat and to tune at will the 
relaxation time of modes of different frequency. This allows one to optimize the time needed to ther- 
malize the system and generate independent configurations. We show how this frequency-dependent 
response can be exploited to control the temperature of Car-Parrinello-like dynamics, keeping at 
low temperature the electronic degrees of freedom, without affecting the adiabatic separation from 
the vibrations of the ions. 
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Solving Hamilton's equations leads to sampling of the 
microcanonical constant-energy distribution, but in real- 
life experiments it is the temperature that is kept con- 
stant. Reproducing this condition in computer simula- 
tions is of great importance for the investigation of a 
large class of physical, chemical and biological problems. 
Several approaches have been proposed to modify Hamil- 
ton's equations in order toperform constant-ternperature 
dynamics (see e.g. Refs. Many of these [l|, 0] rely 

on stochastic methods, which are a natural choice for 
modeling the interactions with an external heat bath, 
and which display excellent ergodic behavior due to their 
random nature. A good thermostat should be able to 
rapidly enforce the correct probability distribution, and 
generate uncorrelated configurations, which arc necessary 
to compute ensemble averages. The efficiency of the ther- 
mostat is particularly important in ab initio simulations, 
because of their high computational cost. The stochastic 
thermostats used so far arc based on Markovian equations 
of motion, and imply no memory of the past trajectory 
of the system. 

Markovian random processes are, however, only a sub- 
set of all possible stochastic processes. Furthermore, the 
Mori-Zwanzig theory ensures that whenever some degree 
of freedom is integrated out, the dynamics of the remain- 
ing degrees of freedom are described by a non-Markovian 
Langevin equation, with a finite-range memory func- 
tion 0-01 • Hence, in the quest for a better thermostat, 
and considering the thermostat as arising from a set of 
bath variables whose effect is integrated out, it is natural 
to explore the effect of using a non-Markovian Langevin 
equation to perform constant-temperature molecular dy- 
namics. In this Letter we will show that, by using colored 
noise, it is possible to influence in a different manner 
the different vibrational modes of the system. There- 
fore the thermostat can be adjusted to the system under 
study, and its performance optimized in a precise and 
predictable fashion. This is, to our knowledge, the first 



time that a colored Langevin equation has been employed 
in atomistic simulations. 

An area which would greatly benefit from an improved, 
tunable thermostat is that of Car-Parrincllo (CP)-like, 
extended Lagrangian schemes The idea behind this 
approach is very general, as it applies to any system 
where the forces are the result of an expensive optimiza- 
tion procedure. This process is circumvented by extend- 
ing the dynamical degrees of freedom (DOF), so as to 
include the parameters to be optimized, and introducing 
an artificial dynamics which allows these extra variables 
to be maintained close to the ground state, by adiabatic 
decoupling from the other degrees of freedom. In the pro- 
totypical example of CP molecular dynamics (CPMD) a 
fictitious mass is assigned to the electronic DOF so that 
they can be evolved at the same time as the ionic DOF. 
If the fictitious mass is small enough, the dynamics of the 
electrons are adiabatically separated from the dynamics 
of the ions. Hence, the electrons are kept close to the 
ground state, while the nuclei are evolved at the correct 
temperature. This same technique can be used in clas- 
sical simulations that use polarizable force fields, where 
the electronic DOF describe the charge polarization of 
the system ■ Similar approaches have also been sug- 

gested in the field of rare-events sampling, to separate the 
oscillations of the microscopic degrees of freedom from 
those of a few selected slow reaction coordinates 11 1. 

Controlling the temperature in these CP-like tech- 
niques requires that one acts separately on the ionic de- 
grees of freedom, which must sample the correct canon- 
ical ensemble, and on the variational parameters, which 
must always remain at low temperature to minimize 
the error in the forces Traditional stochastic ther- 

mostats allow for a highly ergodic sampling of all the de- 
grees of freedom, irrespectively of their frequency. This is 
beneficial for the ionic DOFs but causes the breakdown 
of adiabatic separation. For this reason, deterministic 
thermostats of the Nose-Hoover (NH) type 0] have been 
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adopted. However the original NH thermostat has well- 
known ergodicity problems, and the extension to Nose- 
Hoover chains is normally used 0]. This comes though 
at the price of introducing a large number of parame- 
ters, whose effect on the ions dynamics is not easy to 
predict and control. In the following we show that by us- 
ing correlated noise it is possible to tune the coupling of 
a stochastic thermostat with the various degrees of free- 
dom. This allows one not only to use Langevin dynamics 
in CP-likc methods, but also significantly improves the 
sampling of the target ensemble, because the thermostat 
is tailored to the system under study, in a predictable 
and controlled fashion. 

We consider here a system described by coordinates g^, 
momenta pi and masses m^, interacting via a potential 
U{q), where o is the set of g^'s. The colored Langevin 
equations [1,01 read 



<ii{t) = Pi{t)/m, 



* (1) 



where fi = —dU/dqi are the forces, JC{t) is the mem- 
ory kernel and is a vector of independent Gaus- 
sian noises. In order to set the temperature to a cho- 
sen value T, the noise term C,{t) needs to be related to 
the memory kernel by the fluctuation-dissipation theorem 
{Ut)C3{t'))=6,jm^TlC{t~t'). 

The non-Markovian Eqs. might seem at first too 
complex to be used in practical applications. However, 
for a rather general form of the memory kernel, JC [t) = 
^X]fe '^fe^"*''''^''^'"'''' with 7fc > 0, it is possible to rewrite 
Eq. ((T|) in an equivalent Markovian form by introducing 
a set of auxiliary momenta 13|, [IJ] : 

= {h[q{t)],0, ... ,0)^ - As,(t) + Br7,(i). 

Here Sj = (p^, s^i, . . . , Sij^) is a iV + 1 dimensional vec- 
tor, whose first component is the canonical momentum pi 
associated to the i-th DOF, and "q^ is a vector of Gaussian 
white noises, with {riik{t)rijk'{t')) = - t')5kk'- The 

real-valued matrices A and B determine the dynamics 
of Pi, and can be related to K, (t) by extending the argu- 
ments of Ref. as will be discussed elsewhere. 

In order to illustrate some of the effects of using a col- 
ored noise, we study the simple case in which 
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corresponding to the desired canonical ensemble for q and 
p. The memory kernel and its power spectrum are 
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FIG. 1: (color online) The autocorrelation time of the total 
energy for harmonic oscillators of frequency uj [Cf. Equa- 
tion ((5])] is plotted for different values of the thermostat pa- 
rameters. Dark curves correspond to high friction (7"^ = 20 
fs) whereas light ones correspond to a more gentle thermo- 
stat (7^^ = 1 ps). Dotted lines correspond to white noise 
{tp — 0) and full ones to colored noise with Tp = 2 fs. The 
curves are superimposed on the vibrational density of states 
(DoS) for a polarizable force-field simulation of crystalline cal- 
cite, which was obtained from the Fourier transform of the 
velocity-velocity autocorrelation function. For reference, we 
report the shell vibrational modes as obtained from a run 
where we artificially heated the shells to 300 K. 

respectively. Thus the friction 7 determines the intensity 
of the kernel and Tp the autocorrelation time of the noise. 
For the purpose of this work, one can consider S (lu) to 
be a low-pass filter for the noise, which has the cutoff 
frequency Tp^ . Clearly, when — > the white- noise 
limit is recovered. 

We consider the dynamics of a set of harmonic oscilla- 
tors. In this case Eqs. (0) are fully linear, and the auto- 
correlation time for the total energy of an eigenmode of 
frequency uj can be explicitly evaluated EES: 
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We take th (w) as a measure of the time needed for the 
thermalization of each individual normal mode. For a 
white noise {rp ~ 0), th decreases with to until it reaches 
a plateau at th — I/7, while for Tp ^ 0, the autocorre- 
lation time has a minimum at w = /{2tp) and grows 
quadratically thereafter. By properly adjusting Tp, one 
can select which modes are going to be maximally cou- 
pled with the thermostat, and thus reduce the coupling 
of the thermostat to the fastest modes (see also Fig.[T|). 

We next consider the application of the colored-noise 
thermostat, with the parameters of Eq. ([3]), to classical 
MD simulation using a polarizable force field. Here the 
electronic DOF are represented by charged shells, bound 
with harmonic potentials to the corresponding atomic 
cores. We couple a colored-noise thermostat to the ions, 
at the target temperature, and choose the filtering time 
Tp in such a way that the impact on the electronic DOF is 
minimal. At the same time, we apply a zero temperature 
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FIG. 2: (color online) Shell temperature (Tg) for calcite as 
a function of the thermostat parameters. Simulations have 
been performed for the discrete series of values indicated by 
arrows on the horizontal axes. The points are joined by con- 
tinuous lines, for clarity sake. In both panels we distinguish 
the strength of the ion thermostat by the line color. Darker or 
lighter (red or blue in the onhne version) curves correspond re- 
spectively to a strong (7"^ = 20 fs) or mild (7"^^ ~ 1 ps) fric- 
tion. In panel (a) we plot Ts against tf, and we choose two ex- 
treme values of the shell friction, 7^ ^ = 1 ps and 7g ^ — 50 fs, 
which are represented respectively with full and dashed lines. 
In panel (b) we plot the dependence of Ts versus 73. Here 
full and dashed lines correspond respectively to a physically 
meaningful filter {rp = 2 fs) and to white noise {tf — 0). 

thermostat of friction 7s to the electrons. This latter 
thermostat is memory-less, so that it couples optimally 
with the fast electronic modes. Such a simulation scheme 
amounts to a non-equilibrium dynamics, in which heat 
is injected into the ionic DOF and systematically sub- 
tracted from the electronic ones. In spite of the stochastic 
nature of these equations it is still possible to introduce 
a conserved quantity that measures the accuracy of the 
integration. This can be obtained by accumulating the 
change in kinetic energy due to the thermostat [i , If 



However, at variance with Refs. 0, EB], the conservation 
of this quantity does not rigorously measure the sampling 
accuracy. 

As an example we consider the simulation of crystalline 
calcite, modeled by a polarizable force field [17|]. The 
Ca^"*" ions are treated as non-polarizable, while the polar- 
ization of the COg" anions is described by a charged shell 
attached to each oxygen. The thermostats are applied to 
the non-polarizable ions and, in the case of the oxygens, 
to the center of mass of the system formed by the ion 
plus its shell. Meanwhile, the electronic temperature is 
controlled by the damping of the velocity of the shells 
relative to the partner O ions. The vibrational density of 
states in the absence of any thermostat can be used as an 
approximate guide to the choice of the colored thermo- 
stat parameters (see Fig. [1]) . In real- life, anharmonicity 
will introduce some coupling between the normal modes, 
so that deviations from the predictions of Eq. (O are ex- 
pected. However, at least in the case of quasi-harmonic 
modes, they will most likely reduce th (w). Thus, one can 
safely use the analytical estimate to tune the thermostat 
parameters beforehand, without having to perform time- 
consuming tests on the real system. 

We simulated [l^ a box containing 96 CaCOa units, 



with a timestep of 1 fs, performing NVT runs with target 
temperature T — 300 K. We performed systematic tests 
by varying Tp, 7 and 75 (Fig. |2|. The averages have 
been computed from 1 ns-long runs, where we discarded 
the first 100 ps for equilibration. Within a large range 
of parameters, the procedure performs as expected: the 
temperature of the shells remains below a few K, and 
the ions equilibrate to the desired temperature. As Tp is 
set to a value different from zero, the heat transferred to 
the electronic DOF is reduced. However, some care must 
be taken in choosing the friction 7s, because the shell 
thermostat can induce a small drag on the ions which 
results in an ionic temperature lower than desired, if not 
compensated by a high thermostat strength 7. Since th 
does not decay fast enough to zero for cj > t^^ , one must 
choose a low cutoff frequency in order not to heat up the 
shells. As a consequence, the relaxation time for high- 
frequency phonons increases, making the effects of shell- 
induced drag more pronounced. However, the thermostat 
can be systematically improved by adding more degrees 
of freedom, so as to obtain a more sharply defined filter, 
as we will show below. 

Thermostatting on ab initio CPMD is more challeng- 
ing. Since wavefunctions are not atom-centered, the cou- 
pling of the dynamics of the electronic DOF to the ions 
is stronger than in the shell-model case, and the presence 
of high-frequency components in the noise quickly heats 
up the electrons. Furthermore, because of the expense 
of ab initio CPMD, it is mandatory to have fast equili- 
bration and sampling. We will show that both problems 
can be solved thanks to the tunability and predictability 
of our scheme. As a test example, we ran simulations of 
a single heavy water molecule in vacuum, using a stan- 
dard literature setup (see Fig. |3] and Ref. [lH). We ran 
several independent trajectories for a total of 90 ps, start- 
ing from ionic configurations equilibrated at 300 K and 
from wavefunctions quenched to the Born-Oppenheimcr 
surface We have used Eq. (0) with 5 extended mo- 
menta and fitted A and B in order to obtain a short, 
optimal response time over the ionic degrees of freedom, 
and an abrupt increase in the region corresponding to 
electronic modes [see inset of Fig. ^h)] . We then com- 
pare this case with results from a massive Nose-Hoover- 
chains simulation 0, [23| • In both cases the strength of 
the thermostat is such that the underlying dynamics of 
the ions is severely altered. 

With the present, very conservative choice of parame- 
ters the drift in electronic energy is negligible for both 
thermostats. In Figure [3] we plot the autocorrelation 
function of the squares of the normal modes. The in- 
tegral of these functions measures the time required to 
lose memory of the initial configuration. It is evident 
that the use of an optimized colored-Langevin thermo- 
stat dramatically reduces this time. 

The thermostat we have presented offers a number 
of advantages. It can be used in CP-like, extended- 
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path-integrals MD. We believe that this is only a first 
example, and that colored noise will find many other ap- 
plications in a variety of computational problems. We 
thank Dr. G. Tribello for helping us in the simulation of 
calcium carbonate and for carefully reading the paper. 
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FIG. 3: (color online) Autocorrelation functions for the 
squares of (a) the symmetric stretching and (b) the bend- 
ing modes of a heavy water molecule in vacuum, performed 
in the NVT ensemble at T = 300 K. We use a fictitious mass 
H = 200 a.u., and a time step of 4 a.u., in order to minimize 
the errors on the forces [2y|. The Nose-Hoover thermostat 
with chain length 4 has been used, and its mass chosen so as 
to maximize the coupling to the stretching mode. The NH 
correlation functions (lighter lines, blue in the online version) 
are highly oscillating and decay very slowly. The shading 
highlights the curve's envelope. In contrast using the new 
thermostat (darker lines, red in the online version) we find a 
much sharper decay, which in the case of the stretching re- 
quires an enlarged scale to be appreciated [inset of panel (a)] . 
In the inset of panel (b) we show the relation between th and 
oj for our thermostat. The parameters have been optimized 
to obtain a sharp decay of the response for frequencies above 
the stretching mode. 

Lagrangian simulations, and it is also much faster in 
reaching equilibrium than the Nose-Hoover thermostat. 
This is particularly relevant when performing expensive, 
ab initio simulations, but any problem which requires av- 
eraging over uncorrelated configurations of the system 
can greatly benefit from the enhanced relaxation time. 
The optimal parameters of the simulation can be easily 
estimated before the run is started. Here, in the difficult 
case of a molecule in vacuum, we have been able to re- 
duce the correlation time down to a fraction of a picosec- 
ond. An additional advantage is that the exact propaga- 
tor in the case of zero force is obtained easily [l6|] , which 
makes the implementation simple and robust, at vari- 
ance with Nose-Hoover chains [3| which requires a high 
order integrator to ensure accurate trajectories Fi- 
nally, the introduction of highly tunable, non-Markovian 
thermostats in molecular dynamics simulations lays the 
foundations for the development of optimal sampling al- 
gorithms, which can be of great benefit in free-energy 
techniques, or when one must treat systems with a broad 
vibrational spectrum, which is the case for instance in 
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